REPORT  DOCUMENTATION  PAGE 


Form  Approved 
OMB  No.  0704-0188 


Public  reporting  burden  for  this  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources, 
gathering  and  maintaining  the  data  needed,  and  completing  and  reviewing  this  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this 
collection  of  information,  including  suggestions  for  reducing  this  burden  to  Department  of  Defense,  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and 
Reports  (0704-0188),  1215  Jefferson  Davis  Highway,  Suite  1204,  Arlington,  VA  22202-4302.  Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  law,  no  person 
shall  be  subject  to  any  penalty  for  failing  to  comply  with  a  collection  of  information  if  it  does  not  display  a  currently  valid  OMB  control  number.  PLEASE  DO  NOT  RETURN  YOUR 
FORM  TO  THE  ABOVE  ADDRESS. 


1.  REPORT  DATE  (DD-MM-YYYY)  2.  REPORT  TYPE 

May  2013  T echnical  Paper 


4.  TITLE  AND  SUBTITLE 

A  3D  Unstructured  Mesh  Euler  Solver  Based  on  the  Fourth-Order  CESE 
Method 


3.  DATES  COVERED  (From  -  To) 
May  2013-June  2013 


5a.  CONTRACT  NUMBER 

In-House 


5b.  GRANT  NUMBER 


5c.  PROGRAM  ELEMENT  NUMBER 


6.  AUTHOR(S) 

Bilyeu,  David;  Yu,  John;  Cambier,  Jean-Luc 


5d.  PROJECT  NUMBER 


7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

Air  Force  Research  Laboratory  (AFMC) 

AFRL/RQRS 
1  Ara  Drive. 

Edwards  AFB  CA  93524-7013 


9.  SPONSORING  /  MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 

Air  Force  Research  Laboratory  (AFMC) 

AFRL/RQR 
5  Pollux  Drive 

Edwards  AFB  CA  93524-7048 


12.  DISTRIBUTION  /  AVAILABILITY  STATEMENT 

Distribution  A:  Approved  for  Public  Release;  Distribution  Unlimited.  PA#  13285 


5e.  TASK  NUMBER 


5f.  WORK  UNIT  NUMBER 

Q0A5 


8.  PERFORMING  ORGANIZATION 
REPORT  NO. 


10.  SPONSOR/MONITOR’S  ACRONYM(S) 


11.  SPONSOR/MONITOR’S  REPORT 
NUMBER(S) 

AFRL-RQ-ED-TP-2013-114 


13.  SUPPLEMENTARY  NOTES 

Conference  paper  for  the  AIAA  CFD  conference,  San  Diego,  CA,  June  2013. 


14.  ABSTRACT 

In  this  paper,  the  CESE  method  is  extended  and  employed  to  construct  a  fourth-order,  three-dimensional,  unstructured-mesh 
solver  for  hyperbolic  Partial  Differential  Equations  (PDEs).  This  new  CESE  method  retains  all  favorable  attributes  of  the 
original  second-order  CESE  method,  including:  (i)  flux  conservation  in  space  and  time  without  using  a  one-dimensional 
Riemann  solver,  (ii)  genuinely  multi-dimensional  treatment  without  dimensional  splitting  (iii)  the  CFL  constraint  remains  to  be 
<1,  and  (iv)  the  use  of  a  compact  mesh  stencil  involving  only  the  immediate  neighboring  nodes  surrounding  the  node  where 
the  solution  is  sought.  Two  validation  cases  are  presented.  First  higher  order  convergence  is  demonstrated  by  the  linear 
advection  equation.  Second  supersonic  flow  over  a  spherical  body  is  simulated  to  demonstrates  the  schemes  ability  to 
accurately  resolve  discontinuities. 


16.  SECURITY  CLASSIFICATION  OF: 


a.  REPORT 


b.  ABSTRACT  c.  THIS  PAGE 


Unclassified  Unclassified  Unclassified 


17.  LIMITATION 

OF  ABSTRACT 

18.  NUMBER 

OF  PAGES 

SAR 

16 

19a.  NAME  OF 
RESPONSIBLE  PERSON 

Justin  Koo 


19b.  TELEPHONE  NO 

(include  area  code) 

661-525-5707 

Standard  Form 
298  (Rev.  8-98) 

Prescribed  by  ANSI 


A  3D  Unstructured  Mesh  Euler  Solver  Based  on  the 
Fourth-Order  CESE  Method 


David  L.  Bilyeu  *1,2,  S.-T.  John  Yu  **  1 ,  and  Jean-Luc  Cambier  *2 

1  The  Department  of  Mechanical  and  Aerospace  Engineering,  The  Ohio  State  University, 

Columbus,  OH  43210 

2Air  Force  Research  Laboratory,  Edwards  Air  Force  Base,  CA  93524 

In  this  paper,  the  CESE  method  is  extended  and  employed  to  construct  a  fourth-order,  three-dimensional, 
unstructured-mesh  solver  for  hyperbolic  Partial  Differential  Equations  (PDEs).  This  new  CESE  method 
retains  all  favorable  attributes  of  the  original  second-order  CESE  method,  including:  (i)  flux  conservation  in 
space  and  time  without  using  a  one-dimensional  Riemann  solver,  (ii)  genuinely  multi-dimensional  treatment 
without  dimensional  splitting  (iii)  the  CFL  constraint  remains  to  be  <  1,  and  (iv)  the  use  of  a  compact  mesh 
stencil  involving  only  the  immediate  neighboring  nodes  surrounding  the  node  where  the  solution  is  sought. 
Two  validation  cases  are  presented.  First  higher  order  convergence  is  demonstrated  by  the  linear  advection 
equation.  Second  super  sonic  flow  over  a  spherical  body  is  simulated  to  demonstrates  the  schemes  ability  to 
accurately  resolve  discontinuities. 


I.  Introduction 

Previously  the  CESE  method  has  been  expanded  to  higher  order  for  both  one-dimension,1,  and  two- 
dimensions.  ’  In  this  paper  we  further  extend  the  CESE  method  to  fourth-order  accuracy  in  three-dimensions. 
Similar  to  the  original  second-order  CESE  method,  space  and  time  are  treated  in  a  unified  manner.  Although 
the  present  development  is  fourth-order  the  formulation  to  be  presented  in  the  paper  is  recursive  and  can  be 
straightforwardly  extended  to  a  sixth-,  eighth-order  or  higher  by  including  more  terms  in  the  Taylor  series 
expansion. 

What  makes  the  CESE  scheme  unique  is  due,  in  part,  to  the  unified  treatment  of  space  and  time  as  well  as 
how  the  space-time  domain  is  discretized.  In  general,  the  space-time  domain  is  divided  into  Solution  Elements 
(SEs)  and  Conservation  Elements  (CEs).  In  each  SE  the  primary  unknowns  and  fluxes  are  discretized  and 
represented  by  a  third-order  Taylor  series.  The  CE  area  a  series  of  non-overlapping  regions  that  fill  the  entire 
domain.  These  elements  are  mechanism  in  which  flux  conservation  is  enforced.  The  numerical  integration  is 
aided  by  the  discretized  unknowns  and  fluxes  defined  in  each  SE.  In  general,  CEs  do  not  coincide  with  SEs. 
By  enforcing  flux  conservation  over  each  CE,  a  set  of  algebraic  equations  in  term  of  the  unknown  and  its 
spatial  derivatives  are  derived.  By  solving  the  equations  the  solution  at  the  next  time  step  is  calculated.  As 
will  be  shown  in  the  following  sections,  the  time-marching  calculation  of  the  three-dimensional  fourth-order 
CESE  method  is  explicit. 

The  remainder  of  this  paper  is  organized  as  follows.  Section  II  presents  the  three-dimensional,  fourth- 
order  CESE  method  and  is  divided  into  the  following  sub-sections:  (A)  properties  of  Taylor  series,  (B)  the 
calculation  of  the  fluxes  and  their  derivatives,  (C)  the  calculation  of  the  temporal  derivatives  of  the  conserved 
variables,  (D)  space-time  integration  for  conserved  variables  and  their  even  derivatives  (E)  central  differencing 
approach  for  calculating  the  odd  derivatives  (F)  outline  of  the  numerical  calculation.  Section  III  reports  the 
numerical  results  by  using  the  new  four-order  CESE  method.  At  the  end  of  the  paper,  we  present  the 
concluding  remarks  and  provide  the  list  of  the  cited  references. 
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II.  Numerical  Method 


A  three-dimensional  hyperbolic  equation  is  cast  into  the  following  vector  form: 


dU  dFx 
dt  dx 


dF  y  OF 
dy  + 


dz 


=  0 


(1) 


where  U  =  (u\,u2, . . . ,  umy„  Fx’v’z  =  (fi’v,z,  f2’v’z,  ■  ■  ■  >  /m9’2)*!  and  m  is  the  number  of  equations  in  the 
system.  Each  of  the  equations  in  Eq.  (1)  is  a  divergence-free  condition:  V  ■  h,  =  0,  where  i  =  1, . . . ,  m  and 
h?  =  if?-,  /'(■  ft-  ui)f  is  the  space-time  flux  function.  Here,  the  divergence  operates  in  the  four-dimensional 
Euclidean  space  (a :,y,z,t).  Aided  by  Gauss’  theorem,  the  integral  form  of  Eq.  (1)  becomes 


j)  h i(x,y,z,t)  ■  ds  =  0. 


(2) 


The  CESE  method  is  designed  to  integrate  Eq.  (2)  in  the  space-time  domain.  In  the  numerical  algorithm, 
space  and  time  are  treated  in  a  unified  manner. 

The  geometry  for  a  three-dimensional  CESE  method  is  more  difficult  to  visualize  than  the  one-  and 
two-dimensional  methods  because  the  integration  is  carried  out  in  four- dimensions,  three  spatial  and  one 
temporal.  As  such  the  CE  and  SE  discretization  will  not  be  shown. 

II. A.  The  Taylor  Series  Expansion 

A  Ath-order  Taylor  series  expansion  of  the  conserved  variables  iq  inside  a  SE  can  be  written  in  the  following 
form: 


N-  N- 
N—a  a—b  a—b—c 


N 

Ui(x,y,z,t)  =  ^2  Y  E 

a— 0  6—0  c— 0 


da+b+C+dUi  Axa  A yb  Azc  Atd 


dxadybdzcdtd  al 

d= 0  w 


6! 


dl 


(3) 


where  Ax  =  a :  —  Xj,  Ay  =  y  —  yj,  A z  =  z  —  Zj,  and  At  =  t  —  tn.  The  subscript  j  and  the  superscript  n 
are  respectively  the  spatial  and  temporal  anchoring  point.  For  the  fourth-order  CESE  method  a  third-order 
Taylor  expansion  in  space  and  time  is  employed,  N  =  3.  Since  the  CESE  scheme  does  not  assume  the 
symmetry  property  of  derivatives,  e.g.  gxgy  an<i  dydx  are  ^w0  distinct  variables,  Eq.  (3)  is  not  initially 
compatible  with  this  assumption.  In  order  to  use  Eq.  (3)  within  our  scheme  the  mixed  derivatives  needs  to 
be  averaged.  For  example 


d3 


Ui 


d3 


Ui 


_ .if _ 

dx2dz  3  V  dxdxdz 


d3Ui 

dxdzdx 


d3 


Ui 


-i 


dzdxdx J 

Moreover,  a  derivative  of  Ui  can  also  be  expressed  by  a  Taylor  series  expansion.  All  derivatives  of  iq  in 
Eq.  (3)  can  be  succinctly  expressed  by  the  following  Taylor  series  expansion: 


dcur 


dx1  dyJ  dzK  dtL 


A—  A - 

A  A— a  a—b  a—b—c 

=  EEE  E 

dxI+adyJ+bdzK+cdtL+d  a! 


dBi 


Axa  A yb  Azc  Atd 


bl 


c! 


d\ 


(4) 


where  C  =  I  +  J  +  K  +  L,A  =  N  —  C,  and  B  =  C  +  a  +  b  +  c  +  d.  Obviously,  the  Taylor  series  expansion 
of  m,  Eq.  (3),  is  a  special  case  of  Eq.  (4)  with  A  =  N  and  C  =  0.  Similarly,  the  fluxes,  /f  ’v’z,  and  their 
derivatives  inside  a  SE  are  also  discretized  by  the  Taylor  series  expansion: 


dcf*'y’z 


dx1  dyJ  dzK  dtL 


A-  A- 

A  A— a a—b a—b—c 

=  EEE  E 

a— 0  6—0  c— 0  d— 0 


dBn 


x,y,z 


Axa  A yb  A  zc  A td 


dxI+adyJ+bdzK+cdtL+d  a!  bl  cl  dl 


(5) 


where  C  and  B  have  the  same  definitions  as  that  in  Eq.  (4).  When  C  =  0,  Eq.  (5)  is  the  Taylor  series 
expansion  of  the  flux  f-'v'z  itself. 

The  Taylor  series  coefficients  listed  in  Eqs.  (4)  and  (5)  are  the  unknown  variables  that  need  to  be 
calculated  as  part  of  the  Higher  Order  CESE  method. 

In  the  following  derivation  it  will  be  shown  that  the  only  independent  variables  are  the  conserved  variables 
and  their  spatial  derivatives.  This  requires  that  the  fluxes  and  their  derivatives  as  well  as  the  temporal 
derivatives  of  the  conserved  variables  are  functions  of  the  conserved  variables  and  their  spatial  derivatives. 


DISTRIBUTION  A. 
Approved  for  public  release; 
distribution  unlimited. 


2  of  15 


American  Institute  of  Aeronautics  and  Astronautics 


II. B.  Flux  Terms 


In  a  previous  paper,'  it  was  demonstrated  that  the  flux  terms  and  their  spatial  and  temporal  derivatives 
can  be  expressed  as  functions  of  the  conserved  variables.  This  is  accomplished  through  the  Jacobian  and 
its  derivatives  or  through  a  generalized  Leibniz  rule. 1,1  In  short  both  of  these  methods  demonstrated  this 
by  recognizing  that  the  flux  functions  and  their  derivatives  can  be  written  as  functions  of  the  conserved 
variables,  and  their  derivatives.  The  Jacobian  method  is  more  generic  and  as  such  will  be  used  to  show  this 
relationship. 

Aided  by  the  chain  rule,  the  first  derivatives  of  fz'v'z  can  be  represented  as 


dfj’y’Z  _  y-  df*’v’z  dm 
9’Li  ^  dui  chi') ' 


(6) 


where  ’Ll  =  x,y,  z,  and  t.  On  the  right  hand  side  of  Eq.  (6),  8f^’v’z/8ui  is  the  Jacobian  matrix,  which 
can  be  easily  derived  from  the  relationship  between  the  fluxes  and  the  conserved  quantities.  For  the  second 
derivatives,  we  have 


g2fx,v,z  m  gfx,y,z  g^  m  m  g2 fx,y,z  g^  g^ 

dV1d*2  ~  ^  dui  <9’Lk9’L2  +  ^  ^  8uidup  d’Li  d'K2  ’ 

i=i  i=i  p= l  r 


(7) 


where 


(x,x),  (x,  y),  (x,z),  (x,f), 

\&  )=  ( y ’  A  A  A  A  A A 

(z,x),(z,y),(z,z),(z,t), 
and  (t,t) 

The  second  term  on  the  right  hand  side  of  Eq.  (7),  82 'v'z / du[dup  is  a  m  x  m  x  m,  matrix,  which  can  be 
readily  derived  based  on  the  governing  equations.  For  the  third  derivatives,  we  have 

d3f*'y’z  _^dftv'z  d3Ul 

d^id^2d^3  gUl  d'f?  id^  2dty  3 

<92/f ,v,z  /  d2ui  dup  d2ui  dup  d2ui  dup  \ 

^  ^  duidup  \d^id^2  d^3  +  d’LiS’Ls  d^2  +  d^2d^3  dVi)  +  ^ 

l—l  p= 1  1  x  7 

m  m  m  ,-.3  »x,y,z  000 

®  Jj  ®ui  ®uv  ®uq 
2^  2^  duidupduq  d^i  8^2  8^3  ’ 

(’Ll,'L2,<L3)  = 

(x,  x,  x),  (x,  x,  y),  (x,  x,  z) ,  (x,x,t),  (x,y,x),  (x,y,y),  (, x,y,z ),  (. x,y,t ),  (x,z,x),  ( x,z,y ),  (x,z,z),  ( x,z,t ), 
(y,  x,  x),  (y,  x,  y),  (: y ,  x,  z),  (y,  x,  f),  (y,  y,  x),  (y,  y,  y),  (y,  y,  z),  (y,  y,  t),  (y,  0,  x),  (y,  y),  (y, z,  z),  (y,  z,  f), 

(z,  x,  x),  (z,  x,  y),  (z,  x,  z),  (z,  x,  t),  (z,  y,  x),  (z,  y,  y),  (z,  y,  z),  (z,  y,  f),  (z,  z,  x),  (z,  z,  y),  (z,  z,  z),  (z,  z,  t), 

[x,t,t),  (y,t,t),  ( z,t,t ), 

The  last  term  on  the  right  hand  side  of  Eq.  (8),  83  f*’v,z /8uidupduq  is  a  mxmxmxm  matrix,  which 
can  be  readily  derived  based  in  the  definition  of  f^’v,z  as  functions  of  tq.  Aided  by  Eqs.  (6-8),  we  can  relate 
all  derivatives  of  f:z’v’z  to  the  derivatives  of  iq. 

II. C.  Temporal  Derivatives 

As  with  the  two-dimensional  fourth-order  CESE  scheme  '  the  Cauchy-Kovalewski  procedure  for  non-linear 
equations  is  used  to  relate  the  temporal  derivatives  of  the  conserved  variables  to  the  conserved  variables  and 
its  spatial  derivatives. 
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To  begin  we  let  the  Taylor  series  satisfy  Eq.  (1)  at  point  (j,n) 


( d^y  =  _  (djzy  _  (dfiy  _  (dfiv 

\dt)j  V  dx  J.  \  dy  )  j  V  dz  ). 


(9) 


Each  term  in  the  above  equation  is  a  coefficient  in  a  Taylor  series.  For  conciseness  the  super-  and  sub-scripts, 
i.e. ,  j  and  n,  indicating  the  space-time  location  are  dropped  in  the  following  equations.  Aided  by  Eq.  (9), 
the  second-order  derivatives  of  m  involving  first-order  time  differentiation  are  readily  available  by  assuming 
that  the  following  algebraic  equations  are  valid: 


d2Ui 

d2ff 

d2n 

d2f* 

d2Ui 

d2  fi 

d2ff 

d2f * 

dxdt 

dxdx 

dydx 

dzdx ’ 

dydt 

dxdy 

dydy 

dzdy ’ 

d2Ui 

d2fi 

d2ff 

d2fz 

d2ut 

d2ff 

d2n 

d2f* 

dzdt 

dxdz 

dydz 

dzdz ’ 

dtdt 

dxdt 

dydt 

dzdt1 

Even  though  the  final  relation  in  Eq.  (10)  has  a  temporal  derivative  on  the  right  hand  side  of  the  equation 

Q 2  fx 

it  is  still  a  function  of  spatial  derivatives.  As  an  example  we  look  at  the  first  term  on  the  RHS  dxgt  which 
can  be  expressed  as 


02/f  . .  V  d/f  ff  d2/f  dutdup 

dxdt  ^  dui  dxdt  ^  ^  duidup  dx  dt 

l—l  l  —  l  p=  1  1 

=  _y 'dfl  ( d2f f  d2f f  d2f*\  d2f f  dm  (dJl  dJl  dff\ 

dui  \dxdx  dydx  dzdx)  duidup  dx  \  dx  dy  dz  ) 


This  procedure  can  also  be  applied  to  the  other  terms  in  the  relation  in  Eq.  (10). 

For  the  third-order  derivatives  of  iq  involving  first-order  time  differentiation,  we  assume  that  the  following 
equations  are  valid: 


d3m 

3  ( 

'  Q2fi 

d2/y 

d2f f 

dxd'S’dt 

Kdxdx 

dydx 

dzdx 

d3Ui 

^  ( 

'd2  ft 

d2fi 

d2f t 

dyd^dt 

d<5>  \ 

Kdxdy 

dydy 

dydy 

d3Ui 

9  ( 

f d2f i 

d2  fi 

d2ff 

dzd^dt 

d$  ' 

\dxdz 

dydz 

dzdz 

d3Ui 

d 

(  d2  St 

d2  fi 

d2n 

dtdtdt 

dt 

\dxdt 

dydt 

dzdt 

(11) 


where  'I'  =  (x,  y,  z,  t).  Equation  (11)  has  terms  with  multiple  temporal  derivatives  on  the  RHS  which  can  be 
expressed  in  terms  of  spatial  derivatives  using  the  same  approach  used  for  Eq.  (10). 

Essentially,  Eqs.  (9-11)  assume  that  the  coefficients  of  Taylor  series  expansion  for  iq  and  satisfy 

the  additional  higher-order  equations,  which  are  readily  obtained  by  applying  spatial  differentiation  to  the 
Euler  equations  Eq.  (1).  The  procedure  is  recursive  and  can  be  readily  extended  to  higher-order  derivatives 
of  Ui. 

To  recap,  all  the  first-,  second-,  and  third-order  derivatives  of  iq  and  /;T"/,Z  involving  any  order  of  time 
differentiation  can  always  be  replaced  by  relationship  formulated  in  terms  of  spatial  derivatives  of  uy.  This 
is  achieved  by  the  following  steps:  (i)  the  additional  equations  shown  in  Eqs.  (9-11),  (ii)  Eqs.  (6-8),  and 
(iii)  the  chain  rule  as  shown  in  Eqs.  (6-8).  As  such,  the  independent  variables  in  the  fourth-order  CESE 
method  include  only  the  conserved  variables  Ui  with  i  =  1,2, ...  ,m  and  their  spatial  derivatives.  Table  1 
lists  the  primary  unknowns  of  the  fourth-order  CESE  method. 
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Table  1:  The  unknowns  of  a  three-dimensional,  fourth-order  CESE  method. 


Even-order  variables  Odd-order  variables 


Ui 

{^i)x 

(ui)y 

i)z 

ij^i)xxx 

(v>i)xxy 

i)xxz 

(jJ’i)xy 

i)xyx 

{Ui)Xyy 

i)xyz 

i^i)xz 

i)xzx 

(^i)xzy 

(y,i)xzz 

('ti i)yx 

(V'ijyxx 

( Ui)yxy 

(u i)yxz 

{ui)yy 

{Ui)yyx 

(' ui)yyy 

(Ui)yyZ 

{ui)yz 

i^i)yzx 

{ui)yzy 

(V’i)yzz 

i^i)  zx 

i^i)  zxx 

(ifj)  zxy 

(' Ui)zxz 

( ui)zy 

i^i)  zyx 

(' ui)zyy 

( Ui)Zyz 

( Ui)zz 

fa i)zzx 

(' Ui)zzy 

(' Ui)zzz 

II. D.  Space-Time  Flux  Conservation 

In  this  section,  the  procedure  of  integration  for  space-time  flux  conservation  is  illustrated.  The  space-time 
integration  is  used  to  calculate  the  conserved  variables  and  their  even  derivatives.  Furthermore,  the  method 
presented  in  this  section  is  applicable  to  any  even  derivative,  As  a  reminder  the  integration  is  carried  out  in 
four  dimensions,  three  space  and  one  time.  To  visualize  the  volume  in  which  the  integration  takes  place  it 
is  easiest  to  think  of  it  as  a  cube  where  each  surface  is  really  a  volume.  This  volume  has  sides  and  as  well  as 
a  top  and  bottom  region.  To  more  easily  present  the  integration  it  will  be  split  up  into  two  parts  with  the 
first  being  the  “side”  hyperplanes  and  the  second  being  the  “top  and  bottom”  hyperplanes.  A  hyperplane  is 
a  generic  term  for  any  shape  with  n-1  dimensions  that  can  split  an  n  dimensional  object.  In  two-dimensions 
a  hyperplane  would  be  a  line  and  in  three-dimensions  a  hyperplane  would  be  a  surface. 

The  flux  through  the  side  and  bottom  faces  are  calculated  from  the  known  solution  at  the  previous  time 
step  at  neighboring  points  of  the  current  solution  point.  Therefore,  the  calculation  for  fluxes  through  side 
and  bottom  surfaces  are  explicit.  The  flux  through  the  top  surface  is  a  function  of  the  solution  at  the  new 
time  step. 

Flux  through  Side  Hyperplane 


Figure  1  shows  part  of  the  side  hyperplane.  For  a  tetrahedral  mesh  there  would  be  6  of  these  hyperplanes 
in  one  BCE  giving  a  total  of  24  side  hyperplanes  for  the  entire  CE.  The  points  Pi,  i  =  1,2,3  are  at  the 
current  time  step  while  points  P/,  i  =  1,  2,  3  are  at  the  next  time  step.  Point  Pi  and  P[  have  the  same  spatial 
coordinates  but  different  temporal  location. 


P' 

Pi 


Fig.  1:  The  side  face  associated  with  a  solution  point 
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To  evaluate  the  integral  in  Eq.  (2)  we  calculate  the  normal  via 


i 

j 

k 

t 

H 

bo 

1 

2/2  -  2/1 

N 

to 

1 

b* 

£2  — 

23  -  Xl 

2/3  -  2/i 

z3  -  Zi 

H  — 

x\  —  Xi 

y[  -  yi 

z[  -  Zi 

A- 

By  applying  Eq.  (12)  to  Eq.  (2) 

yields: 

x,y,z  N  N—aN—a—bN—a—b—c 

xr  sr  sr  sr  \' 

dff 

nxi 

/' 

=  [ Nx  i,  Nv  j,Nz  k,0t\ 


(12) 


dxadybdzcdtd  a\b\c\d\  JV4.d 


(x  -  XjT  ( y  y0)b  {z  -  zjf  (t  -  nddVs\  (13) 


where  VgD  represents  the  integration  over  the  side  hyperplane.  This  is  the  integration  will  have  to  be 
repeated  for  each  BCE.  By  noticing  that  the  temporal  integration  is  decoupled  from  the  spatial  integration 
we  can  integrate  the  temporal  terms  separately.  We  can  also  represent  the  volume  through  a  parametric 
equation  with  two  independent  variables,  e.g.  x  =  x\  +  (22  —  x\ )u  +  (23  —  xi)v.  This  leads  to  the  equation, 


N-  N-a- 

x,y,z  N  N—aa—b  b—c 

EEEE  E 

x%  a— 0  b— 0  c— 0  d— 0 


dfl 


N 


(3) 


dxadybdzcdtd  a\b\c\(d  +  1)!  \  2  /  Jo  Jo 


A t_\ 


d-\- 1  pi  pi  — 


[(x2  -  Xi)u+  (X3  —  Xi)  V  +  XI  -  Xj]a 


[(2/2  -yi)u+  ( 2/3  -yi)v  +  yi-  yj]  [{z2  -  Zi)u+  (z3  -  zi)  v  +  zi  -  Zj]c  dvdu, 


where 


i\r(3)  = 


22  -  2l  y 2-  y  1  Z2  -  -21 

X3  -Xl  2/3  -  1/1  Z3  -  Z 1 


(14) 


(15) 


The  integral  in  Eq.  (14)  will  be  integrated  for  select  combinations  of  a,  b  and  c  in  Appendix  A. 
Flux  through  the  Top  and  Bottom  Hyperplane 


To  proceed,  we  present  the  calculation  of  the  flux  through  the  bottom  hyperplane  of  the  CE  centered  at 
point  Pi.  The  bottom  hyperplane  of  the  CE  is  located  in  the  preceding  time  step  n  —  1/2.  The  calculation 
of  the  flux  through  the  bottom  hyperplane  of  the  CE  requires  integration  over  four  triangular  bipyramids. 
A  triangular  bipyramid  is  a  three-dimensional  object  consisting  of  two  tetrahedrons  that  share  a  common 
face.  Figure  2  is  a  diagram  of  a  portion  of  the  top  and  bottom  hyperplane.  This  shape  is  called  a  triangular 
bipyramid  which  is  composed  of  two  tetrahedrons  that  share  a  common  face.  In  this  figure  the  common 
face  is  shown  in  red  while  the  blue  edges  makes  up  the  tetrahedron  associated  with  the  neighboring  cell  and 
the  black  edges  are  associated  with  the  central  cell.  The  point  Pi  is  the  current  solution  point,  point  P4  is 
the  neighboring  cell  center  and  the  point  Pj  is  the  neighboring  solution  point.  When  integrating  over  the 
bottom  hyperplane  the  Taylor  series  is  expanded  from  Pj  but  when  integrating  over  the  top  hyperplane  the 
Taylor  series  is  expanded  from  Pi.  Since  the  geometry  of  the  top  and  bottom  hyperplanes  is  the  same  only 
the  derivation  for  the  bottom  hyperplane  will  be  detailed. 
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Pa 


Fig.  2:  The  bottom  and  top  hyperplane 

To  evaluate  the  integral  in  Eq.  (2)  over  the  bottom  and  top  hyperplane  we  need  the  normal  which  is 


N  = 


i  j  k  t 

x2  -X!  y2-  2/1  £2  -  zi  h~  h 

X3  ~Xi  1/3  -  2/1  Z3  -  zi  t3-  ti 

X4  -  XI  2/4  -  2/1  Zi  -  Z\  ti-  ti 


X2  —Xl  2/2  -  2/1  Z2  -  z  1 

X3  -  Xl  2/3  -  2/1  Z3  —  Z 1 

X4  -  Xl  2/4  -  2/1  Zi  —  Zi 


(16) 


It  should  be  noted  that  the  normals  in  the  i,j  and  k  direction  are  all  zero.  Applying  Eq.  (16)  to  Eq.  (2) 
yields: 


N  N-aN-a-b 

EE  E 

a— 0  b— 0  c=0 


da+b+c  nt 


dxadybdzc  a\b\cl  Jv 


Sri  /  (*  -  ( V  -  yj)b  (z  -  Zj)cdv^D, 

lOlCl 


(17) 


where  n*  is  the  normal  component  in  the  t  direction.  It  should  be  noted  that  there  is  no  time  integration 
in  Eq.  (17)  because  all  points  inside  of  the  volume  VgD  are  at  the  same  time  step.  To  integrate  Eq.  (17) 
we  represent  the  volume  into  its  parametric  form.  After  applying  the  parametric  transformation  Eq  (17) 
becomes 


^^alVg-6  da+b+cUi  |J|  „t  [1 


fl—U  pl—u  —  v 


a— 0  6—0  c— 0 


dxadybdzc  a\b\c\  J0  J0  J0 


x2  —  xi)  u  +  (x3  —  Xi)  v  +  {xi  —  Si)  w  +  xi  —  Xj\a 


[(2/2  -  2/1)  u  +  (2/3  -  2/i)  x  T  (2/4  2/1)  w  +  yi  -  2 lj]b  [(z2  -  zi)u+  ( z3  -  zi)  v  +  (z4  -  zi)  w  +  zi  -  Zjf^dwdvdu 

(18) 

where  |  J\  =  det(iV).  This  leads  to  the  simplihcation  of  \  J\nt  =  N 4.  The  integral  in  Eq.  (18)  will  be  integrated 
for  select  combinations  of  a,  b  and  c  in  Appendix  B. 
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To  progress  an  even  derivative  to  the  next  time  step  we  need  to  apply  either  Eq.  (18)  or  Eq.  (14)  to  each 
hyperplane  in  every  BCE  associated  with  the  current  cell  (j,  n)  and  then  take  the  sum  of  all  integrations. 
For  a  tetrahedral  mesh  each  cell  has  four  BCEs  and  each  BCE  has  three  “side”  hyperplanes,  one  “top” 
hyperplane  and  one  “bottom”  hyperplane. 


II. E.  First  Derivatives 

In  this  section  we  detail  the  procedure  used  to  calculate  the  first  derivatives  of  the  conserved  variables. 
The  method  used  is  very  similar  to  the  one  used  in  the  second-order  scheme,  the  major  difference  is  that  a 
third-order  Taylor  series  is  used  instead  of  a  first-order  Taylor  series.  A  brief  outline  used  to  calculate  the 
first  derivatives  is  detailed  below,  for  a  complete  derivation  please  see  a  second-order  derivation.  ’  As  a 
reminder  the  superscript  and  subscripts  on  the  outside  of  a  bracket  represent  the  location  where  the  Taylor 
series  coefficients  are  calculated  at. 

As  a  preface  we  will  lay  out  the  necessary  geometry.  First  consider  the  current  cell  at  location  jo  and 
the  surrounding  cells  denoted  by  jr  for  r  =  1,2,...,  7Vb ,  where  Nb  is  the  number  of  neighbors.  For  a 
tetrahedral  mesh  Nb  would  be  4.  These  cells  exist  at  both  the  current  time  step,  n,  and  the  previous  time 
step  n  —  1/2  giving  them  the  space-time  location  of  (/r,ro)  r  =  0,  ...,-ZVh  for  the  current  time  step  and 
( jr ,  n  —  1/2)  r  =  0, . . . ,  Nb  for  the  previous  time  step.  There  are  two  important  locations  in  each  cell,  (i)  the 
centroid  denoted  by  j* ,  (ii)  the  Taylor  expansion  point  denoted  by  j+.  The  derivation  below  is  applicable 
to  any  type  of  mesh,  i.e.  tetrahedral,  hexagonal,  prism,  and  four-sided  pyramid,  but  for  simplicity  we  will 
restrict  our  examples  to  tetrahedral  meshes.  In  what  follows  is  an  outline  detailing  each  step  required  to 
calculate  the  first  derivatives. 


1.  Determine  the  location  of  the  Taylor  series  expansion  points  using  either  the  c-r  scheme  detailed  in 
or  the  edge  based  derivative  (EBD)  scheme.'1 

2.  Perform  Taylor  series  expansions  from  (j’o  ,  n)  to  each  (j+,  n)  r  =  1, . . . ,  Nb. 


A  A— a A—a—b 


=  EE  E 


a— 0  6= 0  c— 0 


& 


a-\-b-\-cn 


dxadybdzc 


Jo 


(y-y)  (y-y)  (y~y) 

a\b\c\ 


(19) 


By  applying  Eq  (19)  to  each  of  the  surrounding  cells  we  obtain  Nb  equations  with  three+iVb  unknowns 
per  governing  equation.  The  unknowns  are  the  first  derivatives  of  the  conserved  variables  and  the 
values  of  the  conserved  variables  at  the  expansion  points,  [u*(j+ ,n)]jx  j  =  1, . . . ,  A/. 

3.  Next  the  unknowns  at  [rt*(j+,n)]Jx  are  approximated  by  a  Taylor  series  expansion  from  j*  at  the 
previous  half  time  step  to  j/-  at  the  current  time  step,  i.e.  [u*(j+,  n)]"x  ~  \u'Aj+,  n)]?x  ^  where 

JQ  Jr 


i  n— 1/2 


A  A— a  A—a—b  A—a—b—c 


■V.nJir  =EE  £  E 


a— 0  6=0  c— 0 


d= 0 


dBux 


dxadybdzcdtd  J  ■  x 


(y-y)  (y-y)  (y-y) 

a\b\c\d\ 


t)' 


(20) 


By  substituting  Eq.  (20)  into  Eq.  (19)  and  moving  the  unknowns  to  the  LHS  we  get 


[(u«)*]"x  A^.+  +  [(ui)y}”xAyj+  +  [(«*)z]"xA^.+  = 


Ko+.»)]r,/2-EE  E 


A  A— a A—a—b 


f  da+b+cui  \n  (Ay)  (Ay)  (Ay)  (21) 


a=0  6=0  c= 0 


\  dxadybdzc J  ,x 
Jo 


alblcl 


where  Ax ^  =  x £  -  x. x  ,  A  =  yg  -  ,  A  zg  =  z-+  -  z^  and  (a,  6,  c)  ±  (1,  0, 0),  (0, 1,  0),  (0,  0, 1). 

4.  By  applying  Eq.  (21)  to  neighboring  cells  we  have  now  have  three  unknowns  and  Nb  equations  which 
leads  to  an  overdetermined  system.  This  is  actually  advantageous  because  it  allows  us  to  determine 
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multiple  solutions  and  then  select  the  most  appropriate  solution.  The  possible  solutions  are  determined 
by  grouping  Eqs.  (21)  is  sets  of  three  to  generate  a  potential  solution.  For  example  a  tetrahedral  mesh 
has  four  surrounding  cells  and  would  be  grouped  as  [(1, 2, 3),  (1,  2, 4),  (2,  3, 4)].  Taking  the  first  triplet 
from  the  set  would  result  in  the  following  equation: 


Aaq+  A  y1+  Az1+ 

( Ui )i1) 

Ax2+  Ay2+  A  z2+ 

{Ui)y] 

_Ax3+  Ay  3+  Az3+_ 

{ Ui)(z\ 

K(1+>™)]ix  1/2 

^"~\A —  CL 

—  A^a= 0  2^6=0 

— a — b  j 

2^c= o  1 

(  d^+'ui  ' 

n 

(Axi+y  (Avi+)b  (Azi+y 

^  dxadybdzc ^ 

h? 

n 

a\b\c\ 

K(2+,  n)]£r1/2 

^-~\A  \-^A—a 

Z-/a=0  2-/ 6=0 

\-^A—a—b  j 
2^c=0  ' 

(  da+b+cUi  ' 

(■ Ax2+ y  (■ av2+  ) b  (Az2+  y 

^  dxadybdzc ^ 

io 

n 

a!b!c ! 

K(3+>™)]3^1/2 

^^A — cl 
~  2_/a— 0  2^6=0 

\^A — cl — b  i 

2^c=o  ‘ 

( d°-+b+cUi ' 

(  A:c3+ )“  (  A2t3+ )  6  (Az3+ r 

^  dxadybdzc ^ 

a!b\c\ 

for  (a,  b,  c)?  (0,1,0),  (1,0,0),  (0,0,1). 


(22) 


Using  Eq  (22)  on  all  of  systems  of  equations  will  generate  multiple  solutions  to  (zt,)x ,  ( Ui)v  and  (ui)z. 
The  superscript  (1)  for  the  unknowns  in  Eq  (22)  represents  the  first  possible  solution  for  ( Ui)x ,  {v,i)y 
and  {ui)z. 

5.  The  final  step  is  to  weight  the  possible  solutions  together  to  get  the  optimum  solution.  Any  of  the  pre¬ 
viously  derived  weighting  schemes  are  applicable  to  the  higher  order  CESE  method  with  out  requiring 
any  changes.  In  the  present  investigation  the  W2f'  and  S21  schemes  are  used. 


II. F.  Order  of  Calculation 

In  this  subsection  we  list  the  order  in  which  each  derivatives  is  calculated.  Since  both  the  even  and  odd 
derivatives  requires  information  from  the  current  time  step  the  order  in  which  the  various  derivatives  are 
updated  is  important  to  ensure  that  the  procedure  is  explicit.  We  will  limit  this  procedure  to  a  fourth-order 
scheme  but  the  this  method  can  be  easily  expanded  to  sixth-,  eighth-  or  even  higher  orders  of  accuracy.  For 
a  fourth-order  scheme  this  is  the  order  in  which  the  values  are  calculated. 

1.  Temporal  derivatives  of  the  conserved  variables  and  the  fluxes  and  their  spatial  and  temporal  deriva¬ 
tives. 

2.  Second  derivatives,  e.g.  (ui)xx,(ui)xy,(ui)xz 

3.  Third  derivatives,  e.g.  (ut)xxx:  ('Uji)xxyi  (.'Uijxxz 

4.  Conserved  variables,  Ui 

5.  First  derivatives,  ( Ui)x,(ui)y,(ui)z . 


Second  Derivatives 


To  proceed,  we  apply  all  possible  second  derivatives  with  repetition  to  Eq.  (1)  to  yield  the  following  additional 
equations: 

<92h, 

V  '  d^id'^2  =  0l  =  {(x>x)i(xiV)>(x>z)’  (y,x),(y,y),(y,z),  (z,  x),  (z,  y),  (z,  z)}  . 

Recall  that  h  =  (/* ,  ff,  ff ,  ui)  is  the  space-time  flux  vector.  Aided  by  the  Gauss  theorem  in  the  three- 
dimensional  space-time  domain,  the  above  differential  equations  are  recast  into  the  following  integral  equa¬ 
tions: 


d'S’idi&2 


■  ds  =  0  (^r  ,^2)  =  {{x,x),(x,y),(x,z),  (y,  x),  (y,  y),  (y,  z),  (z,x),(z,y),(z,z)}  .  (23) 


9  of  15 


DISTRIBUTION  A. 
Approved  for  public  release; 
distribution  unlimited. 


American  Institute  of  Aeronautics  and  Astronautics 


In  the  setting  of  the  fourth-order  CESE  method,  ( Ui)xx  is  represented  by  using  the  first-order  Taylor  series 
expansion  with  (uj)txx,  ( Ui)xxx ,  and  ( Ui)yxx  as  the  coefficients.  The  integration  can  be  straightforwardly 
performed  based  on  the  original  second-order  CESE  method.  Essentially,  we  invoke  the  original  CESE  nine 
times  to  calculate  the  second  order  spatial  derivatives  at  the  solution  point  (j.  n) . 

Third  Derivatives 


Once  the  second-derivatives  are  calculated,  the  third-derivatives  of  it,  listed  in  Table  1  are  calculated  by 
central-differencing  the  second  derivatives  of  14.  By  slightly  modifying  the  procedure  outline  in  Section  II. E 
we  obtain  the  third  derivatives.  For  example  to  determine  ( Ui)xyx ,  ( Ui)xyy ,  ( Ui)xyz  we  would  obtain 


Ax1+ 

Ayi+ 

A  zi+ 

~(u  A1) " 

)xyx 

Ax2+ 

Ay2+ 

A  z2+ 

( Ui)xyy 

_Ax3+ 

Ay3+ 

A  z3+_ 

{^i)^cyz 

'[K)*y(l +,n)T1x1/2 

^~^A  ^^A — a, 
~  2^a= 0  2^6=0 

Y^A — cl — 6 

2^c=o 

f  da+b+cm 

r 

(Aa:l+)a(A!/l+),,(Ai:l+)C" 

K dxadybdzc 

r 

a!6!c! 

[«M2+,n)£-1/2 

^-~\A  \-^A—a 

2-/a=0  6—0 

\-^A—a—b 

2^c=0 

(  da+b+cu.i 

(Ax2+y  (AV2+)b  (Az2+T 

^  dxadybdzc 

'jo* 

r 

a!6!c! 

[(U'M3+,  n)]”"1/2 

Y^A  ^"~rA  —  CL 

2_/a— 0  2-/ 6=0 

\  ~\  A — a — 6 
2^c=o 

(  aa+b+cui 

(Ax3+Y  (AV3+)b  (Az3+T 

^  dxadybdzc 

iox 

a!6!c! 

for  (a,  b,  c)±  (0,1,0),  (1,0,0),  (0,0,1). 


Which  simplifies  to 


Ax\+ 

Ax2+ 

Ax3+ 


Ay  1+ 
Ay2+ 
Ay3+ 


Azi+ 

{.'U"i)xyx 

A  z2+ 

( Ui)xyy 

A  z3+_ 

(W’ij^xyz 

M*„(l +,n)]n1~1/2 
[K)*„(2  +:n)]n2~1/2 
[KM3+,  n)]”-1/2 


[(ui)xy]™x 

[{ui)xy]™x 

[(ui)xy]™x 

Jo 


(24) 


Eq.  (24)  is  nearly  identical  to  the  equation  used  to  determine  the  first  derivatives  in  the  second-order  scheme. 
This  procedure  is  then  used  to  determine  all  other  third  derivatives. 


Conserved  Variables 


The  next  step  is  to  update  the  conserved  variables,  iq.  For  this  step  you  use  the  method  outlined  in  Section 
II. D.  The  only  modification  required  is  to  set  N  =  3  in  Eqs.  (14)  and  (18).  Although  the  integration  over 
the  top  surface  uses  solutions  from  the  current  time  step  the  scheme  is  still  explicit.  This  is  possible  because 
the  second  and  third  derivatives  are  all  ready  known  and  the  integration  of  the  first  derivatives  is  zero.  This 
leaves  iq  as  the  only  unknown  and  is  easily  separated.  The  integration  of  the  first  derivatives  are  zero,  not 
the  first  derivatives  them  self. 


First  Derivatives 

The  first  derivatives  are  calculated  using  the  procedure  listed  in  Section  II. E. 


III.  Results  and  Discussion 

To  assess  the  accuracy  of  the  three-dimensional,  unstructured  hyperbolic  solver,  we  consider  the  following 
benchmark  problems:  (i)  the  advection  equation  to  assess  the  convergence  rate,  and  (ii)  the  Euler  equation 
to  simulate  supersonic  flow  over  a  spherical  blunt  body. 
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III. A.  Advection  Equation  for  Convergence 

The  advection  equation,  Eq.  (25),  is  used  to  verify  the  order  convergence  of  the  new  scheme. 


du  du  du  du 

~7T7  +  o-x  q — b  ay  — — b  az  —  =  0 
at  ox  oy  dz 

A  solution  is  assumed  to  take  the  form 

u  =  sin(oxa;  +  avy  +  azz  +  att), 


(25) 


(26) 


the  domain  is  a  cube  of  length  2-7t  and  periodic  boundary  conditions  are  imposed.  Under  these  assumptions 
the  wave  speeds  are  ax  =  ay  =  az  =  1  and  at  =  —  (X  +  ay  +  a%).  The  simulation  is  allowed  to  run  for 
a  non-dimensional  time  of  25  which  allows  the  wave  to  progress  through  the  domain  25/(27 r)  times.  To 
determine  the  rate  of  convergence  the  L2  norm  is  calculated  and  compared  against  a  characteristic  length. 
The  characteristic  length  is  taken  to  be  (Volume /ncell)-1^3,  where  ncell  is  the  number  of  cells.  The  result 
of  the  simulation  is  shown  in  Figure  3.  A  best  fit  line  was  taken  and  the  convergence  rate  was  determined 
to  be  4.3  and  2.6  for  the  fourth  and  second-order  schemes  respectively. 


4th 


2  nd 


X 


Fig.  3:  Convergence  test  for  the  3D  convection  equation 


III.B.  Supersonic  Flow  over  a  Blunt  Body 

To  determine  the  schemes  capability  to  resolve  discontinuities  supersonic  flow  over  a  sphere  is  considered. 
The  properties  of  interest  are  the  post  shock  density,  shock  standoff  distance  and  the  shock  profile.  For  the 
simulation  the  working  fluid  is  air  with  specific  heat  ratio  1.4  and  a  gas  constant  287.15  J/kg  K.  The  free 
stream  Mach  number,  pressure  and  density  are  3.0,  1  bar,  and  1.23  kg/m3  respectively.  The  radius  of  the 
sphere  is  0.5  meters.  The  simulation  was  run  for  5  thousand  iterations  at  an  average  CFL  number  of  0.533. 
The  post  shock  density  was  calculated  to  be  4.613  kg/m3  which  compares  favorably  to  the  analytical  value 
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of  4.744  kg/m3.  The  shock  standoff  distance  was  approximated  to  be  0.11  which  also  compares  favorably  to 
the  value  of  0.10  predicted  by  Ambrosio  and  Wortman’s11  relation, 

A  =  0.1431?  e3'24/M°°,  (27) 


where  R  is  the  radius  of  the  sphere  and  is  the  free  stream  Mach  number. 

The  profile  of  the  shock  around  the  sphere  is  obtained  by  the  relationship  generated  by  Billig  and  is 
equal  to 


x  =  R  +  A  —  Rc  cot2  6 


(  y2 tan2  9 


(28) 


where  A  is  the  shock  stand  off  distance  given  by  Eq.  (27),  9  is  the  Mach  angle  is  equal  to  sin  1(1/M00)  and 
Rc  is  the  radius  of  curvature  and  is  equal  to 

Rc  =  1. 14377  e(0-54/(M°o-i)1'2_ 

A  numerical  Schlieren  image  is  shown  in  Figure  4.  This  figure  shows  the  shock  profile  generated  by  the 
simulation  and  the  points  represent  shock  profile  generated  by  Eq.  (28).  The  shock  profile  predicted  by 
Eq.  (28)  agrees  favorably  with  the  numerical  results. 

IV.  Concluding  Remarks 

This  paper  details  the  derivation  and  validation  of  a  new  fourth-order  CESE  method  used  to  solve  three- 
dimensional  hyperbolic  PDEs  on  unstructured  tetrahedral  meshes.  The  formulation  was  found  to  retain  all 
favorable  features  of  the  original  second-order  CESE  method,  including  (i)  the  use  of  the  most  compact  mesh 
stencil  involving  only  the  immediate  neighboring  nodes  of  the  central  node  where  the  unknowns  are  sought, 
(ii)  the  stability  constraint  of  the  four-order  CESE  method  remains  to  be  CFL  <  1,  and  (iii)  completely 
explicit  operation  in  the  time  marching  calculation.  To  demonstrate  the  capabilities  of  the  new  method, 
two  test  cases  were  reported:  (i)  a  sinusoidal  wave  modeled  by  the  advection  equation,  used  to  confirm 
higher-order  convergence  (ii)  Mach  3  air  flow  around  a  spherical  blunt  body  used  to  determine  the  schemes 
ability  to  accurately  resolve  discontinuities.  Future  work  will  investigate  the  effects  of  various  re-weighting 
schemes  on  the  accuracy  of  the  results  and  the  development  of  a  three-dimensional  Navier-Stokes  solver. 

A.  Side  Integration 

This  section  will  list  selected  integrations  for  the  integral  in  Eq.  (14).  To  start  we  will  define  A x;  =  Xi—Xj , 
Ay.;  =  y;  —  yj,  A z;  =  z;  —  Zj  for  *  =  1,2,  3.  With  this  definition  Eq.  (14)  becomes 


r» 

^a^bjC 


1 —u  3 


0  JO 


n  [a^ 


ATgU  ■ 


A4>{  (1 


3= 1 


dvdu, 


where  T  =  {a:,  y,  z}  and  ui  =  {a,  6,  c}. 

When  a  =  1  and  b  =  c  =  0  the  integration  is 


C\, 


0,0 


1  3 

=  sX>*- 


i= 1 


When  a  =  b  =  1  and  c  =  0 


rs 

■‘-1,1,0 


1 

24 


T:  (AxiAyi)  +  ^  Axi  ^2  AVi 


2=1 


2=1 


2=1 
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^1 


Fig.  4:  Convergence  test  for  the  3D  convection  equation 
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When  a  =  b  =  c  =  1 


1  3  3  3  3 

£i,i,i  =  120  [  5Z  Axi  A^  y  +  2  y  ( Axi  Ay,  Az,)  + 

2—1  2=1  2=1  2=1 

3  3  3  3  3  3 

y~l  ( AxiAyi )  y  Azj  +  y  (AxiAzi)  y  Ay,  +  y  {AytAz,)  y  Az,;  . 

2=1  2  =  1  2=1  2=1  2=1  2=1 

Given  these  three  formulas  it  is  possible  to  easily  derive  all  other  integrals  for  a  fourth-order  scheme.  For 
example  to  obtain  the  formulation  for  £2  o  o  one  may  substitute  Axi  for  Ay,  in  the  formulation  for  C\  1  0. 

B.  Bottom  and  Top  Integration 

This  section  will  list  selected  integrations  for  the  integral  in  Eq.  (18).  To  start  we  will  define  Ax,  =  Xi~Xj , 
Ayi  =  yi  —  yj ,  Azj  =  z,  —  Zj  for  i  =  1,  2, 3.  With  this  definition  Eq.  (18)  becomes 


p  1  p  1  — 14  pl  —  U  —  V  3 

/  /  /  +  ASfr^v  +  AfyJ4w  +  (1  —  u  —  v  —  w) 

Jo  Jo  Jo  j=1  *- 


pbt 

n.  h  n 


dwdvdu, 


where  T  =  {x,  y ,  z }  and  u  =  { a ,  6,  c}. 

When  a  =  1  and  6  =  c  =  0  the  integration  is 


nbt  _  1 

^1,0,0  -  24 


y  Ax*. 

i=l 


When  a  =  6  =  1  and  c  =  0 

'4  4  4' 

y;  (A xlAyi)  +  y  Aay  y  Ay i 
.2=1  2=1  2=1 

When  a  =  b  =  c  =  1 


M,1,0  — 


120 


^  4  4  4  4 

1,1  =  [  y  y  Ay,  y  Az.;  +  2  y  (AxiAyzAzi)  + 

2=1  2=1  2=1  2=1 

4  4  4  4  4  4 

y  ( AxiAyi )  y  Az, + y  (a^az,)  y  a^  +  y  (AayAzd  y  az, 

2=1  2=1  2=1  2=1  2=1  2=1 


As  with  the  derivation  in  Appendix  A  the  other  integrals  required  for  a  fourth-order  scheme  are  easily 
obtained  through  substituting.  For  example  Zl|f  j  2  is  obtained  by  substituting  Az^  for  A  ay  in  the  formulation 

of 
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